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The field of frustrated magnetism has been enriched significantly by the discovery of various 
kagome lattice compounds. These materials exhibit a great variety of macroscopic behaviours rang¬ 
ing from magnetic orders to quantum spin liqnids. Using large-scale exact diagonalization, we con¬ 
struct the phase diagram of the S' = 1/2 J 1 -J 2 kagome Heisenberg model with z-axis Dzyaloshinskii- 
Moriya interaction D^. We show that this model can systematically account for many of the experi¬ 
mentally observed phases. Small J 2 and Dz can stabilize respectively a gapped and a gapless spin liq¬ 
uid. When J 2 or Dz is substantial, the ground state develops a Q = 0, 120°antiferromagnetic order. 

The critical strengths for inducing magnetic transition are Df ~ 0.1 Ji at J 2 = 0, and Jf ~ 0.4 Ji 
at Dz = 0. The previously reported values of Dz and J 2 for herbertsmithite [ZnCu 3 (OH) 6 Cl 2 ] place 
the compound in close proximity to a quantum critical point. 

PACS numbers: 75.10.Jm, 75.10.Kt, 75.30.Kz 


In frustrated magnetism m. the kagome lattice has 
become the paradigmatic system of choice for studying 
novel spin-liquid phases that result from geometric frus¬ 
tration and quantum fluctuation [1]. For example, the 
S' = 1/2 antiferromagnetic (AFM) kagome Heisenberg 
model with only the nearest-neighbour (NN) exchange Ji 
is magnetically disordered with a spin correlation length 
less than one lattice spacing [SHU]- The ground state 
is characterized as a gapped topological spin liquid with 
a finite triplet gap 0.05 — 0.15 Ji [TUkTT] and a nonzero 
topological entanglement entropy da eg. Advances in 
computational methods and theoretical ideas have led to 
a deeper understanding and classification of these exotic 
states of matter, which do not break any symmetry and 
sustain emergent fractional excitations. However, a wide 
gulf still separates theory and experiments. 

The celebrated herbertsmithite, ZnCu 3 (OH) 6 Cl 2 [El 
[50], has nearly perfect kagome planes consisting of 
S'=l/2, Cu^+ atoms. Despite a predominant Ji ~ 200 
K, this compound does not develop any long-range mag¬ 
netic order down to T = 50 mK [5TI - f55] . agreeing with 
the NN AFM kagome Heisenberg model. But there are 
disagreements: Without finding any apparent gap down 
to 0.1 meV, neutron scattering suggests gapless exci¬ 
tations [531 Hg. Magnetic susceptibility Xm shows an 
upturn at low temperature mi ng, which is also un¬ 
expected; a gapped spin liquid would otherwise show a 
vanishingly small Xm close to T = 0. To account for 
these discrepancies, Dzyaloshinskii-Moriya (DM) inter¬ 
action [551150] , exchange anisotropy [31] |35] , or quenched 
site dilution [571I301I31I31] have been investigated for 

herbertsmithite. On the other hand, not all kagome com¬ 
pounds are spin liquids. In materials such as Cu(l,3- 
bde) [Cu-benzenedicarboxylate] [SSI |3g, vesignieite 
[BaCu 3 (V 04 ) 2 ( 0 H) 2 ] [57H50] . and Cs 2 Cu 3 SnFi 2 [511 - 


l43] . the ground states is a Q = 0, 120°AFM order, 
which highlights potential interactions beyond the NN 
exchange. Determining the importance of additional cou¬ 
plings is thereby key to a comprehensive understanding 
of the diversified properties in kagome materials [44] . 

In this Letter, we study the effects of the z-axis DM 
interaction and second NN exchange coupling J 2 , 
which are arguably the two most relevant perturbations 
for isotropic kagome compounds. For the first time, we 
construct the phase diagram of the Ji-J 2 -Dz model us¬ 
ing large-scale exact diagonalization (ED) with cluster 
sizes up to = 42 sites. By also examining the excita¬ 
tion gaps and static structure factors, we show that the 
model can sustain various phases including long-range 
AFM order and quantum spin liquids with or without a 
finite spin gap. The distinct ground states among differ¬ 
ent kagome compounds can be systematically accounted 
for with varying interaction strengths of the systems. In 
particular, the reported values of and J 2 for herbert¬ 
smithite indicate a ground state closely proximal to a 
magnetic quantum critical point, where small extra per¬ 
turbations suffice to suppress its long-range magnetism. 
Based on a numerically unbiased method, our study pro¬ 
vides a direct road map for gauging Dz and J 2 in S' = 1 /2 
materials with nearly prefect isotropic kagome structures. 

Model and Method - We consider the J 1 -J 2 kagome 
Heisenberg model with DM interactions: 

H = Ji^ S, • S, + J 2 ^ S, • S, + ^ D,, • (S, X S,), 
<ij> <ij> 

( 1 ) 

where the first two terms are respectively superexchange 
interactions between NN and second NN S = 1/2 sites. 
The third DM-interaction term originates from relativis¬ 
tic spin-orbit coupling and is nonzero when lattice in- 
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version symmetry is absent [451146) . Here we focus on 
the z-axis component of the DM interaction D^- = 
using the convention that Dz > 0 when all links i —)■ j 
are oriented clockwise [see inset of Fig. 1(a)] [37]. We 
neglect the in-plane component Z?||, as it is reported to 
be smaller than Dz in materials of interests and also re¬ 
ducible to second order in Dy/Ji with a spin basis ro¬ 
tation |29l HSl [49]. Throughout the paper we consider 
antiferromagnetic couplings (Ji, J 2 > 0 ) and set Ji = 1. 

We solve Eq. 0 systematically by numerical diag- 
onalization on clusters of size iV up to = 42. The 
Hamiltonian matrix is constructed utilizing translational 
symmetry and conservation (SD) . The re¬ 

sulting sparse matrix eigenvalue problem is solved by the 
Krylov-Schur algorithm as implemented in SLEPc m 
and PETSc |S3| libraries. The cluster choices and calcu¬ 
lation details are given in the Supplemental Material. 

Phase Diagram - We first establish the phase diagram 
of the Jx-J^-Dz model. Without Dz and J 2 , the system is 
magnetically disordered. A finite Dz or J 2 could support 
a Q = 0, 120°AFM ground state with spins lying in the 
cry-plane [see inset of Fig. 1(a)]. Therefore, we proceed to 
map out the phase boundary between the quantum AFM 
state and the magnetically disordered region by studying 
the transverse spin-spin correlation function |29j : 

94 

S:UQ) ( 2 ) 

1.1 

Here, R/,j are unit-cell positions and a, b are intra-unit- 
cell site indices. <S'®^(Q) represents the elements of a 
3x3 matrix and peaks at Q = 0. The largest matrix 
eigenvalue = S 1200 (N) corresponds to the 120° AFM spin 
arrangements, and its classical value is equal to 1 with 
the pre-factor choice of Eq. ([^ [33|. Spontaneous spin 
symmetry breaking can be identified on finite-size clus¬ 
ters by a linear 1/VN extrapolation (SO] |53l El] . When 
the extrapolated Si 2 o° = 'S'i 20 °(-^ = > 0, long-range 

magnetic order develops. 

Figure 1(a) shows the phase diagram obtained by a 
grid interpolation of 11 x 7 = 77 points on the {Dz, J2) 
plane. The generic features are computed with cluster 
sizes TV = 15 — 33, and further N = 36, 39 calculations 
are performed to more precisely locate the phase bound¬ 
ary. The blue region of Fig. 1(a) represents a magnet¬ 
ically disordered ground state {Sf 2 o° < 0)j and the red 
represents the Q = 0, 120°AFM phase {Sf 2 o° > 0)- Fig¬ 
ures 1(b) and 1(c) show linear I/'/N extrapolations of 
S'f^Qo {N) close to the critical transition points along the 
Dz- and J2-8o<.es, respectively. The critical strengths are 
found to be ~ 0.1 at J 2 = 0, and J| ~ 0.4 at Dz = 0. 
We note that Fig. 1(c) shows a more apparent finite-size 
effect with an even-A/odd-iV alternation. With addi¬ 
tional A = 42 calculations in the zero momentum sector, 
the critical strengths when extrapolated independently 
are J|(even-A) ~ 0.32 and J 2 (odd-A) ^ 0.44 at Dz = 0. 
J| is further reduced when Dz is finite, and vice versa. 



FIG. 1. (a) Phase diagram of the Ji-J 2 -Dz model. The 

false-color intensity represents S 1200 . The system is mag¬ 
netically disordered in the blue region, and it develops the 
Q = 0, 120° AFM order in the red. The geometric symbols 
denote estimated values of Dz in various kagome compounds: 
O for herbertsmithite, A for Cu(l,3-bdc), □ for vesignieite, O 
for barlowite, and O for Cs 2 Cu 3 SnFi 2 . (b)-(c) Linear 1/y/N 
extrapolations of Si 2 oo{N) at various Dz and J 2 . 


The Ji-Dz and J1-J2 models have been separately in¬ 
vestigated before. In particular, 11° = 0.1 at J 2 = 0 
was also reported by previous ED studies (211 jSOj- Ji 
at Dz = 0 has been computed by a number of meth¬ 
ods [55H55] . with the reported critical strength rang¬ 
ing from 0.2 to 0.7. Classically, a positive infinitesimal 
J 2 would favour a Q = 0, 120°AFM long-range order 
(J| = 0+) [nniinil- Our finding of J| ~ 0.4 highlights the 
role of quantum fluctuation in destabilizing magnetism. 
Nonetheless, when J 2 becomes substantial, the quantum 
Q = 0, 120° AFM ground state can be stabilized. 

Materials Relevance - We next discuss the relevance of 
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our phase diagram to different S = 1/2, Cu-based ma¬ 
terials with (nearly) perfect isotropic kagome structures, 
as denoted by the geometric symbols in Fig. 1(a). In 
Cu(l,3-bdc), the material develops the Q = 0, 120°AFM 
order with a critical transition temperature ~ 2 
K [33135]. Its interaction strengths J 2 ) ~ (0.1, 0.05) 
estimated by first-principles calculations |53| indeed cor¬ 
respond to a positive Sf^o^ ~ 0.05 in our phase dia¬ 
gram. The compound vesignieite also develops the Q = 0 
magnetic order at Tjv ^ 9 K [37HiD] . Its experimen¬ 
tally estimated DM interaction ~ 0.15 Ji renders a pos¬ 
itive Sf 2 oo ~ 0.1. Similarly, Cs 2 Cu 3 SnFi 2 is ordered at 
Ttv ^ 20 K [ 4 III 43 ] . where the reported could be as 
large as 0.25 Ji, leading to ~ 0.15. Interestingly, a 
higher T;v in experiments seems to be correlated with a 
larger positive Sf 2 o° in our phase diagram. 

Herbertsmithite, however, is not magnetically ordered. 
This could mean that its interactions (D^, J 2 ) are below 
the critical strengths (11°, J|) ^ (0.1, 0.4). On the other 
hand, although J 2 is likely small compared to J 2 , is 
reported to be comparable to 11° in this compound. For 
example, electron spin resonance experiment estimates a 
DM-interaction strength ^ 0.08 Ji |3H]- Theoretical fit of 
Xm indicates a Dz ~ 0.1 Ji m- In fact, first-principles 
calculations find {Dz,J 2 ) ~ (0.1,0.05) for both Cu(l,3- 
bdc) and herbertsmithite |62j , whereas the former is mag¬ 
netically ordered but the latter is not. When the system 
resides in close proximity to the boundary of quantum 
phase transition, additional perturbations such as spin- 
space exchange anisotropy mm I63ll66| or quenched 
dilution of Cu sites [33 1301 1331 I33| in herbertsmithite 

could more easily suppress its long-range magnetism. 

More recently, a structurally perfect kagome com¬ 
pound —the barlowite [Cu 4 (OH) 6 FBr]— has been syn¬ 
thesized [67] . This material develops long-range magnetic 
order at T/v = 15 K with a weak ferromagnetic moment. 
In terms of our phase diagram, such a relatively high Tn 
would imply a strong Dz 0.2 Ji, which in conjunction 
with a small Dy can lead to a Q = 0, canted AFM state 
with a net ferromagnetic moment pointing outward from 
the kagome plane [681 ISO]. Future single-crystal mea¬ 
surements can further clarify the nature of barlowite’s 
low-temperature magnetic structure. 

Gapped versus Gapless Spin Liquids - We next address 
the issue of spin gap by focusing on the = 1 

excitation: At = = 1) - = 0)- Here, 

i?o(5'/Qtj,i) is the lowest energy in a given sector, 

and At corresponds to the triplet excitation when spin 
SU{2) symmetry is present. We will consider only even- 
N clusters based on the gap definition. 

We first note that extrapolating the gap on finite-size 
clusters is more difficult. In the disordered region, the 
scaling form is a priori unknown. The result can depend 
largely on the extrapolation function [HIM 133, as well 
as the cluster size and shape 113 . In the ordered phase, 
the energy splitting between quasi-degenerate ground 



FIG. 2. The spin gap At computed on clusters of size N 
along (a) the J 2 -axis {Dz = 0) and (b) the D^-axis (J 2 = 0). 
The corresponding ground state energy per site Eq/N is given 
in (c) and (d), respectively. 

states in quantum antiferromagnets would scale to zero 
as 1/N [TO] 113 , and the linear-dispersing Goldstone 
modes for spontaneous broken symmetries would scale 
as 1/VN [33] [73] [73]. But there may be no clear sepa¬ 
ration between these states in relatively small systems. 
Despite these difficulties, however, the variation of the 
gap with parameters would be robust and distinguish¬ 
able in our numerically exact data. Therefore, instead of 
making precise quantitative statements, we would mainly 
focus on the trends. 

Figure 2 shows the spin gap At and the ground state 
energy per site Eq/N on different size clusters. As shown 
in Fig. 2(a), At is quickly enhanced by J 2 [73], but 
the rate of increase becomes smaller above J 2 0.1 and 
tends to level off with increasing N. If a simple 1/A 
scaling is employed in the disordered region as in pre¬ 
vious ED study m, At would reach its maximum at 
J 2 ~ 0.1 and then decrease monotonically above it. This 
extrapolated gap behaviour agrees well with recent den¬ 
sity matrix renormalization group calculations [59] . im¬ 
plying that the gapped spin liquid phase is most stable 
around J 2 = 0.1. At large J 2 , although the raw data of 
At appear to grow with J 2 across the magnetic phase 
boundary [J 2 (even-A) ^ 0.32)], we note that the abso¬ 
lute value of the ground state energy |Ao| is also increas¬ 
ing [Fig. 2(c)]. The ratio At/|Ao 1 becomes nearly flat 
and decreases systematically with increasing N in the or¬ 
dered regime. The gap would eventually scale to zero in 
the thermodynamic limit, corresponding to spontaneous 
spin symmetry breaking of the ordered ground state. 
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The effect of on the spin gap is quite differ¬ 
ent. As shown in Fig. 2(b), At is rapidly reduced 
in the presence of a small m- In fact, a simple 
1/A extrapolation would indicate that At is already 
zero in the disordered region before reaching the criti¬ 
cal point D'^ ~ 0.1. This shows the possibility of a gap¬ 
less spin liquid ground state induced by spin exchange 
anisotropy [761177]. At Dz > 0.1, the ratio At/|Ao| stays 
flat and again decreases systematically with increasing 

N. The thermodynamic-limit At remains gapless in the 
magnetic phase. 

The above results show that small J 2 and Dz could 
stabilize respectively a gapped and a gapless spin liq¬ 
uid. Based on the magnetically disordered ground state 
and gapless excitations found in herbertsmithite, our 
study suggests in this material a DM-interaction strength 

O. 05 < Dz < 0.1, closely proximal to the quantum critical 
point 11°. With the prevalent observation of gapless exci¬ 
tations in putative spin-liquid phases, our study also im¬ 
plies that the DM interaction is in general non-negligible 
in isotropic kagome systems. 

Static Structure Factor - To connect with neutron scat¬ 
tering experiments, we study the static structure factor: 

5“^(Q) = (3) 

d 

where a, P G {x,y,z}. Neutron scattering also provides 
energy-resolved spectra by measuring the dynamic struc¬ 
ture factor S°‘^{Q,uj), where S'“^(Q) = f dujS°‘^(Q,uj). 
A spin liquid phase would produce continuous or diffu¬ 
sive scattering spectra, whereas an ordered magnet would 
generate sharp, discrete Bragg peaks. 

Figurej^shows the transverse component S^^(Q) com¬ 
puted on the N = 36 cluster. When Dz and J 2 are both 
zero [Fig. |^a)], the spectrum is close to being uniformly 
distributed along the extended Brillouin zone (BZ); the 
first BZ contains little spectral weight m- This suggests 
that spin correlations are predominantly antiferromag¬ 
netic, while correlation lengths are on the order of lattice 
spacing. Due to spin SU{2) symmetry at Dz = 0, the 
longitudinal component 5'^^(Q) is identical to S'“^(Q), 
and both components are zero at the F point. 

When Dz = 0.125, the Q = 0, 120°AFM ground state 
manifests a structure factor that peaks at the midpoints 
of the extended BZ edges [Fig. |^b)]. In this case, S'°°(Q) 
is much weaker than S^^{Q) aX Dz yf 0, and spins mainly 
lie in the xy-plane. When the system is ordered, 
at large Dz and at large J 2 are in general similar, except 
that (i) S^^{Cl = F) is no longer zero in the former, and 
(ii) = K') is further suppressed in the latter. In 

addition, the overall spectra do not undergo a sharp tran¬ 
sition across the critical point or J|. These features 
can be seen in Fig. [^c) that shows high-symmetry line 
cuts of at various J 2 with Dz = 0. 

In herbertsmithite, neutron scattering signals are dif¬ 
fused for all the measured energies between 0.25 to 11 



FIG. 3. Transverse static structure factors S'^^(Q) of the 
Ji-J 2 -Dz model for J 2 = 0 at (a) Dz = 0 and (b) Dz = 0.125. 
The first and extended Brillouin zones are indicated by the 
dashed and solid hexagons, respectively, (c) line cuts 

along high-symmetry points of the Brillouin zone at various 
J 2 with Dz = 0. 


meV (where Ji ~ 17 meV) [IS]. The spectral weight 
is concentrated in the extend BZ but does not peak at 
any specific Q point, although at 0.75 meV additional 
peak appears at the midpoints of the extended BZ edges. 
These results agree with our 5'^“(Q) calculation for the 
magnetically disordered state [Fig. [^a)]. The experi¬ 
mental spectra also contain a small but finite weight at 
the F point, which could result from the DM interaction. 
We note, however, that the experimental intensity inte¬ 
grated up to 11 meV contains only ^ 20% of the total 
spectral weight. A more detailed theory-experiment com¬ 
parison would require direct calculations of S°‘^ (Q,w). 

In conclusion, we have studied the interplay be¬ 
tween further-neighbour exchange and Dzyaloshinskii- 
Moriya interaction on the kagome lattice. The phase 
diagram of the Ji-J 2 -Dz model is shown to contain 
various novel states of matter, including a Q = 0, 
120 °antiferromagnetic long-range order, as well as 
gapped and gapless quantum spin liquids. A small varia¬ 
tion of the parameters near the phase transition bound¬ 
ary could potentially account for the distinct properties 
observed in different kagome materials. The phase dia¬ 
gram thereby serves as a benchmark for determining the 
importance of these additional perturbations. Studying 
the dynamical properties in different parts of the phase 
diagram and making further connection to inelastic neu¬ 
tron or x-ray scattering measurements would be impor¬ 
tant for future research. 
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In this Supplemental Material, we describe the details 
of exact diagonalization for solving the J 1 -J 2 kagome 
Heisenberg model with z-axis Dzyaloshinskii-Moriya in¬ 
teraction. The calculations are performed with full pe¬ 
riodic boundary conditions on clusters of size iV = 3 x 
Ni X N 2 , ranging from TV = 15 to 42. The employed 
cluster geometries are shown in Fig. SI, while the lattice 
basis vectors remain in agreements with Ref. [1]. Where 
possible, clusters with a higher point group symmetry are 
chosen; for example, both N= 27 and 36 clusters have 
the full spatial symmetry of the kagome lattice. Also, 
the preference is given to cluster geometry with the low¬ 
est possible ground-state energy. 

One exception to our single choice of cluster for a given 
N occurs at = 18. This is because the cluster “18a” 
has a higher point-group symmetry, while “18b” has a 
lower ground-state energy. Due to finite-size effects, the 
resulting spin-gap values seem to vary more substantially 
on the iV = 18 clusters, although the effects of J 2 and 
on the gap trends are similar. Therefore, we have con¬ 
sidered both “18a” and “18b”, by showing the averaged 
spin gaps and ground-state energies in the main text. 

After deciding on the cluster size and geometry, we 
utilize translational symmetry and = "Yhi con¬ 

servation to construct the Hamiltonian matrix on a basis 
of translation operator eigenstates with a fixed [2]. 

The ground state resides in the •S'^otai = sector with 
“momentum quantum number” Q = ^7^)^ where 

ki € {0,1,..., A—1}. We need to diagonalize the matrix 
in each momentum sector to obtain the true ground state. 
For N = 42, we restrict ourself to the Q = (0,0) sector, 
which is indeed the ground-state sector close to the mag¬ 
netic transition points for all smaller clusters under study. 
In this case, C' 2 -rotation and Z 2 -spin-inversion symme¬ 
tries are applied to further reduce the Hilbert space size 
from the original size ~538 billion to the manageable 
size of 9,611,848,877 in the C' 2 -rotation even and Z 2 -spin- 
inversion odd sector for diagonalization. 

In our large-scale sparse-matrix diagonalization, we 
use distributed memory parallelization based on message 
passing interfere, as well as the SLEPc [3] library for 
matrix-eigenvalue problem and the PETSc [4] library for 
distributed vector and matrix operations. The Krylov- 
Schur [5] iterative algorithm chosen for the matrix diago¬ 
nalization is an improved variation of the Arnoldi method 
with optimized implicit restarting and ability to obtain 
multiple degenerate eigenstates. 



FIG. SI. Kagome clusters of different sizes employed in the 
exact diagonalization calculations. The black dots denote the 
choice of spin sites. 
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